home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 007 / threed.arc / THREEDS.F77 < prev   
Encoding:
Text File  |  1986-03-28  |  12.3 KB  |  502 lines
  1. C* */ /* THREEDS.F77   15-Feb-86   G. Doughty                     */
  2. C*****************************************************************************/
  3. C                                         */
  4. C Subroutines to plot function of two variables given as a 100 by 100         */
  5. C array of real values.  Current version drives a Televideo 921 terminal     */
  6. C with the graphics option.                             */
  7. C                                         */
  8. C This code was derived from a program I obtained from a graduate chemistry  */
  9. C student at the University of North Carolina at Chapel Hill in 1973.         */
  10. C (I wish I could give his full name but I only remember that his last name  */
  11. C was Hoskins and that he was working for Dr. Peterson at the time.)  I've   */
  12. C only modified the original code as necessary for compilation with         */
  13. C Digital Research's FORTRAN compiler (as well as indenting it to make it    */
  14. C easier to read).  I've also added routines to replace the original plot    */
  15. C functions (PLOT, PICSIZ, ORIGIN) available on the university's computation */
  16. C center computer.  Finally, I've added comments where I've been able to     */
  17. C fathom the operation of the routines (The original program is largely      */
  18. C devoid of program... but then many of the programs I wrote at UNC are      */
  19. C similarly unadorned).                                 */
  20. C                                         */
  21. C*****************************************************************************/
  22.  
  23.     SUBROUTINE THREED(A, N, M, K)
  24. C*****************************************************************************/
  25. C                                         */
  26. C A is a 100 by 100 real array containing the function values at uniform     */
  27. C grid points.  This array need not be completely filled, however.  N and    */
  28. C M define the limits of the first and second (respectively) array index.    */
  29. C Thus the function occupies A(1 - N, 1 - M).  K is a mystery as the momemt. */
  30. C The example main-line that drives this always passes the value of 3 for K. */
  31. C Unfortunately the other arguments are passed in COMMON.  They are ANGA,    */
  32. C ANGB, and HV.  ANGB appears to tilt the 3-D plot and ANGA appears to ro-   */
  33. C tate it (I'm not sure about the order of these operations on the image).   */
  34. C HV seems to be a height parameter which I have been setting to 5.   I      */
  35. C thinks it controls how much perspective affect one gets.  Beware that      */
  36. C greater or lesser values may cause the image not to fit within 640 x 240.  */
  37. C The program works best if the real values are scalled to the range from    */
  38. C -4.0 to +4.0.                                     */
  39. C                                         */
  40. C*****************************************************************************/
  41.  
  42.     COMMON ANGA, ANGB, HV, D
  43.     COMMON DL, DM, DN, CX, CY, CZ, QX, QY, QZ, DD
  44.     DIMENSION H(10), V(10), X(2), Y(2), Z(2), XP(8), A(100,100)
  45.     CALL    PICSIZ(15.,10.)
  46.     DL = -COS(ANGA) * COS(ANGB)
  47.     DM = -SIN(ANGA) * COS(ANGB)
  48.     DN = -SIN(ANGB)
  49.     IF(ABS(DN).NE.1.0) GO TO 10
  50.         WRITE(6,20)
  51. 20        FORMAT(1X,'******* You are attempting to look straight ',
  52. +        'down (or up) at the surface')
  53.         GO TO 2150
  54. C
  55. 10    CONTINUE
  56.     DD = 1.0/SQRT(1.0 - DN * DN)
  57.     X(1) = 0
  58.     X(2) = N
  59.     Y(1) = 0
  60.     Y(2) = M
  61.     T = N
  62.     IF(N .LT. M) T = M
  63.     D = M * M + N * N + T * T
  64.     D = SQRT(D)
  65.     SCALE = 10.0 * D
  66.     CX = -DL * SCALE
  67.     CY = -DM * SCALE
  68.     CZ = -DN * SCALE
  69.     QX = CX + D * DL
  70.     QY = CY + D * DM
  71.     QZ = CZ + D * DN
  72. 2060    CALL MINMAX(A, N, M, Z)
  73.     DOLR = T / (Z(2) - Z(1))
  74.     DO 30 I = 1, N
  75.         DO 30 J = 1, M
  76.         A(I,J) = (A(I,J) - Z(1)) * DOLR
  77. 30    CONTINUE
  78.     Z(1) = 0.0
  79.     Z(2) = T
  80. 2080    CALL CUBE(X, Y, Z, XP, H, V)
  81.     DO 2130 I = 1, 8
  82.         H(I) = ((XP(I) - QX) * DM - (H(I) - QY)*DL) * DD
  83.         V(I) = (V(I) - QZ) * DD
  84. 2100    CALL MINMAX(H, 8, 1, H(9))
  85. 2120    CALL MINMAX(V, 8, 1, V(9))
  86. 2130    CONTINUE
  87.     IF(ANGB .GT. 0.0) GOTO 2140
  88.         T = V(9)
  89.         V(9) = V(10)
  90.         V(10) = T
  91. 2140    CALL ORTHO (X, Y, A, N, M, H, V, K)
  92. C    CALL PLOT(HV + 2.0, 0.0, -3)
  93.     XXXX = HV + 2.
  94.     CALL ORIGIN(XXXX, 0., 1)
  95. 2150    CONTINUE
  96.     END
  97. C
  98.     SUBROUTINE MINMAX(A, N, M, Z)
  99.     DIMENSION Z(1), A(100, 1)
  100. C
  101. 1050    U = A(1, 1)
  102. 1060    V = U
  103. C
  104. 1080    DO 1190 J = 1, M
  105. 1100        DO 1180 I = 1, N
  106. 1120        IF(U - A(I,J)) 1150, 1150, 1130
  107. 1130            U = A(I, J)
  108. 1150        IF(A(I, J) - V) 1180, 1180, 1160
  109. 1160            V = A(I, J)
  110. 1180        CONTINUE
  111. 1190    CONTINUE
  112. 1210    Z(1) = U
  113. 1220    Z(2) = V
  114. 1230    END
  115. C
  116.     SUBROUTINE CUBE(X, Y, Z, XP, H, V)
  117.     DIMENSION X(1), Y(1), Z(1), H(1), V(1), XP(1)
  118. C
  119. 50    L = 0
  120. 70    DO 180 I = 1, 2
  121. 90        DO 170 J = 1, 2
  122. 110        DO 160 K = 1, 2
  123. 130            L = L + 1
  124. 140            CALL ROTATE(X(I), Y(J), Z(K), XP(L), H(L), V(L))
  125. 160        CONTINUE
  126. 170        CONTINUE
  127. 180    CONTINUE
  128.     END
  129. C
  130.     SUBROUTINE ROTATE(X, Y, Z, XP, YP, ZP)
  131.     COMMON ANGA, ANGB, HV, D
  132.     COMMON DL, DM, DN, CX, CY, CZ, QX, QY, QZ, DD
  133.     DK = D/((X - CX) * DL + (Y - CY) * DM + (Z - CZ) * DN)
  134.     XP = CX + DK * (X - CX)
  135.     YP = CY + DK * (Y - CY)
  136.     ZP = CZ + DK * (Z - CZ)
  137.     END
  138. C
  139.     SUBROUTINE ORTHO(X, Y, A, N, M, H, V, K)
  140.     COMMON ANGA, ANGB, HV, D
  141.     COMMON DL, DM, DN, CX, CY, CZ, QX, QY, QZ, DD
  142. C
  143.     DIMENSION X(1), Y(1), H(1), V(1), A(100,100)
  144.     INTEGER UP, DOWN, PEN, P, Q, P1, PO
  145. C
  146. C
  147. C
  148.     END = 1.0 / 16.0
  149. C
  150. C Can use 1/32 or 1/64 for finer interpolation
  151. C
  152.     UP = 1
  153.     DOWN = 0
  154.     SH = HV / (H(10) - H(9))
  155.     SV = HV / (V(10) - V(9))
  156.     MM = M
  157.     NN = N
  158. 80    IF(K - 1) 100, 120, 100
  159. 100        IF(K - 3) 1110, 120, 1110
  160. 120        L = 0
  161.         LD = 1
  162.         DDD = 0.5 * LD
  163. 140        DO 1060 J = 1, M
  164.             Q = 0
  165.             YJ = J
  166. 160            DO 1030 I = 1, NN
  167.             L = L + LD
  168.             XI = L
  169.             CALL CANSEE(A, XI, YJ, N, M, P)
  170.             PEN = UP
  171.             IF(P) 510, 520, 530
  172. 510            CONTINUE
  173. C            IF(Q) 540, 550, 170
  174.             IF(Q) 540, 550, 540
  175. 520            CONTINUE
  176.             IF(Q) 610, 1020, 610
  177. 530            CONTINUE
  178. C            IF(Q) 170, 550, 540
  179.             IF(Q) 540, 550, 540
  180. 540            CONTINUE
  181.             PEN = DOWN
  182.             GO TO 170
  183. 550            CONTINUE
  184.             IF(I .EQ. 1) GO TO 170
  185.             DI = DDD
  186.             TO = L - LD
  187.             T = TO + DI
  188.             P1 = Q
  189. 560            IF(ABS(DI) .LT. END) GO TO 570
  190.                 CALL CANSEE(A, T, YJ, N, M, PO)
  191.                 DI = DI * 0.5
  192.                 IF(PO .EQ. 0) GO TO 565
  193.                 TO = T
  194.                 P1 = PO
  195.                 T = T - DI
  196.                 GO TO 560
  197. 565                T = T + DI
  198.                 GOTO 560
  199. 570            CONTINUE
  200.             T = TO
  201.             IF(P1 * P) 170, 170, 580
  202. 580            CONTINUE
  203. 590            CONTINUE
  204.             ZP = A(L-LD,J) + (T-L+LD)*(A(L,J)-A(L-LD,J))/LD
  205.             CALL ROTATE(T,YJ,ZP,XP,HH,VV)
  206.             HH = ((XP - QX)*DM - (HH - QY) * DL) * DD
  207.             VV = (VV - QZ) * DD
  208.             HH = (HH - H(9)) * SH
  209.             VV = (VV - V(9)) * SV
  210.             CALL PLOT(HH, VV, PEN)
  211. C600            PEN = 5 - PEN
  212. 600            PEN = 1 - PEN
  213.             GOTO 170
  214. C
  215. 610            CONTINUE
  216.             PEN = DOWN
  217.             DI = DDD
  218.             TO = L - LD
  219.             T = TO + DI
  220.             P1 = Q
  221. 620            IF(ABS(DI) .LT. END) GO TO 630
  222.                 CALL CANSEE(A, T, YJ, N, M, PO)
  223.                 DI = DI * 0.5
  224.                 IF(PO .EQ. 0) GO TO 625
  225.                 TO = T
  226.                 P1 = PO
  227.                 T = T + DI
  228.                 GO TO 620
  229. 625                T = T - DI
  230.                 GOTO 620
  231. 630            CONTINUE
  232.             T = TO
  233.             IF(P1 * Q) 600, 600, 590
  234. 170            CALL ROTATE(XI, YJ, A(L, J), XP, HH, VV)
  235.             VV = (VV - QZ) * DD
  236.             HH = ((XP-QX)*DM - (HH - QY)*DL) * DD
  237. 190            HH = (HH - H(9)) * SH
  238.             VV = (VV - V(9)) * SV
  239.             CALL PLOT(HH, VV, PEN)
  240. 1020            Q = P
  241. 1030            CONTINUE
  242. C
  243. C
  244.             L = L + LD
  245.             LD = -LD
  246.             DDD = -DDD
  247. C
  248. 1060        CONTINUE
  249. C
  250. C
  251. 1090        IF(K - 3) 2060, 1110, 2060
  252. C
  253. 1110        CONTINUE
  254. C
  255.         L = 0
  256.         LD = 1
  257.         DDD = 0.5 * LD
  258. 1140        DO 2040 I = 1, N
  259.             XI = I
  260.             Q = 0
  261. 1160            DO 2020 J = 1, MM
  262.             L = L + LD
  263.             YJ = L
  264.             CALL CANSEE(A, XI, YJ, N, M, P)
  265.             PEN = UP
  266.             IF(P) 1510, 1520, 1530
  267. 1510            CONTINUE
  268. C            IF(Q) 1540, 1550, 1170
  269.             IF(Q) 1540, 1550, 1540
  270. 1520            CONTINUE
  271.             IF(Q) 1610, 2010, 1610
  272. 1530            CONTINUE
  273. C            IF(Q) 1170, 1550, 1540
  274.             IF(Q) 1540, 1550, 1540
  275. 1540            CONTINUE
  276.             PEN = DOWN
  277.             GO TO 1170
  278. 1550            CONTINUE
  279.             IF(J .EQ. 1) GO TO 1170
  280.             DI = DDD
  281.             TO = L - LD
  282.             T = TO + DI
  283.             P1 = Q
  284. 1560            IF(ABS(DI) .LT. END) GO TO 1570
  285.             CALL CANSEE(A, XI, T, N, M, PO)
  286.             DI = DI * 0.5
  287.             IF(PO .EQ. 0) GO TO 1565
  288.             TO = T
  289.             P1 = PO
  290.             T = T - DI
  291.             GO TO 1560
  292. 1565            T = T + DI
  293.             GO TO 1560
  294. 1570            CONTINUE
  295.             T = TO
  296.             IF(P1 * P) 1170, 1170, 1580
  297. 1580            CONTINUE
  298. 1590            CONTINUE
  299.             ZP = A(I,L-LD) +(T-L+LD)*(A(I,L) -A(I,L-LD))/LD
  300.             CALL ROTATE(XI, T, ZP, XP, HH, VV)
  301.             HH = ((XP - QX) * DM - (HH - QY) * DL) * DD
  302.             VV = (VV - QZ) * DD
  303.             HH = (HH - H(9)) * SH
  304.             VV = (VV - V(9)) * SV
  305.             CALL PLOT(HH, VV, PEN)
  306. C1600            PEN = 5 - PEN
  307. 1600            PEN = 1 - PEN
  308.             GO TO 1170
  309. 1610            CONTINUE
  310.             PEN = DOWN
  311.             DI = DDD
  312.             TO = L - LD
  313.             T = TO + DI
  314.             P1 = Q
  315. 1620            IF(ABS(DI) .LT. END) GO TO 1630
  316.                 CALL CANSEE(A, XI, T, N, M, PO)
  317.                 DI = DI * 0.5
  318.                 IF(PO .EQ. 0) GO TO 1625
  319.                 TO = T
  320.                 P1 = PO
  321.                 T = T + DI
  322.                 GO TO 1620
  323. 1625                T = T - DI
  324.                 GO TO 1620
  325. 1630            CONTINUE
  326.             T = TO
  327.             IF(P1 * Q) 1600, 1600, 1590
  328. 1170            CALL ROTATE(XI, YJ, A(I,L), XP, HH, VV)
  329.             HH = ((XP-QX)*DM - (HH-QY)*DL) * DD
  330.             VV = (VV - QZ) * DD
  331. 1180            HH = (HH - H(9)) * SH
  332. 1190            VV = (VV - V(9)) * SV
  333.             CALL PLOT(HH, VV, PEN)
  334. 2010            Q = P
  335. 2020            CONTINUE
  336.             L = L + LD
  337.             LD = -LD
  338.             DDD = -DDD
  339. 2040        CONTINUE
  340. 2060        CONTINUE
  341. 2130    CONTINUE
  342.     END
  343. C
  344.     SUBROUTINE CANSEE(Z, XI, YJ, M, N, P)
  345.     COMMON ANGA, ANGB, HV, D
  346.     COMMON DL, DM, DN, CX, CY, CZ, QX, QY, QZ, DD
  347.     INTEGER CUM, CNT, P, POLD
  348.     REAL I, J, II, JJ
  349.     DIMENSION Z(100, 100)
  350.     IR = XI
  351.     JC = YJ
  352.     ZB = Z(IR, JC)
  353.     IF(XI .EQ. IR) GO TO 2
  354.         ZB = Z(IR,JC) + (XI-IR) * (Z(IR+1,JC) - Z(IR,JC))
  355.     GO TO 4
  356. 2    IF(YJ .EQ. JC) GO TO 4
  357.         ZB = Z(IR,JC) + (YJ-JC) * (Z(IR,JC+1) - Z(IR,JC))
  358. 4    CONTINUE
  359.     XEND = 0.0
  360.     DX = 0.0
  361.     YMULT = 0.0
  362.     ZMULT = 0.0
  363.     IF(XI .EQ. CX) GO TO 10
  364.         YMULT = (YJ - CY) / (XI - CX)
  365.         ZMULT = (ZB - CZ) / (XI - CX)
  366.         DX = 1.0
  367.         XEND = M + 1
  368.         IF(XI .LT. CX) GO TO 10
  369.         DX = -1.0
  370.         XEND = 0.0
  371. 10    CONTINUE
  372.     YEND = 0.0
  373.     DY = 0.0
  374.     XMULT = 0.0
  375.     IF(YJ .EQ. CY) GO TO 20
  376.         XMULT = (XI - CX) / (YJ - CY)
  377.         IF(ZMULT .EQ. 0.0) ZMULT = (ZB - CZ) / (YJ - CY)
  378.         DY = 1.0
  379.         YEND = N + 1
  380.         IF(YJ .LT. CY) GO TO 20
  381.         DY = -1.0
  382.         YEND = 0.0
  383. 20    CONTINUE
  384.     CUM = 0
  385.     CNT = 0
  386.     P = 0
  387.     XB = XI
  388.     YB = YJ
  389. 30    CONTINUE
  390.     II = AINT(XB)
  391.     JJ = AINT(YB)
  392.     XSTEP = DX
  393.     YSTEP = DY
  394.     IF(XB .EQ. II) GO TO 40
  395.         IF(DX. LT. 0.0) XSTEP = 0.0
  396.         GOTO 45
  397. 40    IF(YB .EQ. JJ) GO TO 45
  398.         IF(DY .LT. 0.0) YSTEP = 0.0
  399. 45    CONTINUE
  400.     I = II + XSTEP
  401.     J = JJ + YSTEP
  402.     IF(I .EQ. XEND) GO TO 81
  403.     IF(J .EQ. YEND) GO TO 81
  404.         XB = CX + XMULT * (J - CY)
  405.         YB = CY + YMULT * (I - CX)
  406.         IF(DX .LT. 0.0) GO TO 55
  407.         IF(XB .LT. I) GO TO 60
  408. 50            XB = I
  409.             GO TO 65
  410. 55        IF(XB .LT. I) GO TO 50
  411. 60        YB = J
  412. 65        CONTINUE
  413.         ZB = CZ + ZMULT * (XB - CX)
  414.         IR = I
  415.         JC = J
  416.         IF(YB .NE. J) GO TO 70
  417.         IDX = I - DX
  418.         ZS = Z(IR, JC) - DX*(XB-I)*(Z(IDX,JC) - Z(IR,JC))
  419.         GO TO 75
  420. 70        JDY = J - DY
  421.         ZS = Z(IR,JC) - DY*(YB-J)*(Z(IR,JDY) - Z(IR,JC))
  422. 75        CONTINUE
  423.         SGN = 1
  424.         IF(ZB .LT. ZS) SGN = -1
  425.         CUM = CUM + SGN
  426.         CNT = CNT + 1
  427.         IF(IABS(CUM) .EQ. CNT) GO TO 30
  428.         GO TO 90
  429. 81    CONTINUE
  430.     P = 1
  431.     IF(CUM) 84, 86, 90
  432. 84        P = -1
  433.     GO TO 90
  434. 86    CONTINUE
  435.         IF(ZB .LE. CZ) GO TO 90
  436.         P = -1
  437. 90    CONTINUE
  438.     END
  439. C
  440. C DUMMY ROUTINES TO REPLACE PLOT FUNCTIONS
  441. C
  442.     SUBROUTINE PICSIZ(X, Y)
  443.     REAL X, Y
  444.     CHARACTER GFPL*8
  445.     DATA GFPL/' * mcgW1'/
  446.     IF(X.EQ.0.0) GOTO 3001
  447.     IF(Y.NE.0.0) GOTO 3005
  448. 3001    WRITE(6, 3002)
  449. 3002    FORMAT(1X,'D')
  450.     GOTO 3011
  451. C
  452. 3005    WRITE(6,5000)X,Y
  453. 5000    FORMAT(1X,'PICSIZ: X DIMENSION = ',F10.5,'  Y DIMENSION = ',
  454. +    F10.5)
  455.     GFPL(1:1) = CHAR(27)
  456.     GFPL(3:3) = CHAR(27)
  457.     WRITE(6, 3010)GFPL
  458. 3010    FORMAT(1X,A)
  459. 3011    CONTINUE
  460.     END
  461. C
  462.     SUBROUTINE PLOT(X, Y, PEN)
  463.     COMMON ANGA, ANGB, HV, D
  464.     COMMON DL, DM, DN, CX, CY, CZ, QX, QY, QZ, DD
  465.     REAL X, Y
  466.     REAL YSCALED
  467.     INTEGER PEN
  468.     INTEGER IX, IY, IXO, IYO
  469. C    WRITE(6,5001)PEN, X, Y
  470.     IX = IFIX(X * 639./HV)
  471.     YSCALED = 50.*Y
  472.     IF(YSCALED.LT.0.0) YSCALED = 0.0
  473.     IF(YSCALED.GT.239.) YSCALED = 239.
  474.     IY = IFIX(YSCALED)
  475.     IF(PEN.LT.0) GOTO 3111
  476.     IF(PEN.GT.1) GOTO 3111
  477.     IF(PEN.NE.0) GOTO 3110
  478.         WRITE(6, 3100) IXO, IYO, IX, IY
  479. 3100        FORMAT(1X,'M',I3.3,',',I3.3,';L',I3.3,',',I3.3,';')
  480. 3110    IXO = IX
  481.     IYO = IY
  482.     RETURN
  483. C
  484. 3111    CONTINUE
  485.     WRITE(6,5005)
  486. 5005    FORMAT(1X,'D')
  487.     WRITE(6,5001)PEN, X, Y
  488. 5001    FORMAT(1X,'PLOT : Pen = ',I4,'  at location (',F10.5,F10.5,')')
  489.     WRITE(6,5004)CHAR(27)
  490. 5004    FORMAT(1X,A1,'m')
  491.     END
  492. C
  493.     SUBROUTINE ORIGIN(X,Y,PEN)
  494.     REAL X, Y
  495.     INTEGER PEN
  496. C    WRITE(6,5003)
  497. C    WRITE(6,5002)PEN, X, Y
  498. C5002    FORMAT(1X,'ORIGIN : 3rd arg = ',I4,' 1st arg = ',F10.5,
  499. C+    ' 2nd arg = ',F10.5)
  500. C5003    FORMAT(1X,'D')
  501.     END
  502.